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Abstract. Discontinuous Galerkin methods are developed for solving the Vlasov-Maxwell system, methods that are 
designed to be systematically as accurate as one wants with provable conservation of mass and possibly total energy. Such 
properties in general are hard to achieve within other numerical method frameworks for simulating the Vlasov-Maxwell system. 
The proposed scheme employs discontinuous Galerkin discretizations for both the Vlasov and the Maxwell equations, resulting 
in a consistent description of the distribution function and electromagnetic fields. It is proven, up to some boundary effects, that 
charge is conserved and the total energy can be preserved with suitable choices of the numerical flux for the Maxwell equations 
and the underlying approximation spaces. Error estimates are established for several flux choices. The scheme is tested on the 
streaming Weibel instability: the order of accuracy and conservation properties of the proposed method are verified. 
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1. Introduction. In this paper, we consider the Vlasov-MaxweU (VM) system, the most important 
equation for the modehng of coUisionless magnetized plasmas. In particular, we study the evolution of a 
single species of nonrelativistic electrons under the self-consistent electromagnetic field while the ions are 
treated as uniform fixed background. Under the scaling of the characteristic time by the inverse of the plasma 
frequency u}~^, length by the Debye length Ad, and electric and magnetic fields by —mcujp/e (with m the 
electron mass, c the speed of light, and e the electron charge), the dimensionless form of the VM system is 



where the equations are defined on Q ~ fl^ x fl^, where x d denotes position in physical space, and 
^ G in velocity space. Here /(x,<^,t) > is the distribution function of electrons at position x with 
velocity ^ at time t, E(x, i) is the electric field, B(x, is the magnetic field, p{x,t) is the electron charge 
density, and 3(x,t) is the current density. The charge density of background ions is denoted by pi, which is 
chosen to satisfy total charge neutrality, Jj^ (p(x, — pi) dx — 0. Periodic boundary conditions in x-space 
are assumed and the initial conditions are denoted by /o = /(x,^, 0), Eq — E(x,0) and Bq — B(x, 0). We 
also assume that the initial distribution function /o(x, ^) e n Ll{n^), i.e., is in a Sobolev space of 

order m and is integrable with finite energy in ^-space, where LP^{fl^) = : |'0I^(1 + < oo}. 

The initial fields Eo(x) and Bo(x) are also assumed to be in H"^{ilx). 

The VM system has wide importance in plasma physics for describing space and laboratory plasmas, 
with application to fusion devices, high-power microwave generators, and large scale particle accelerators. 
The computation of the initial boundary value problem associated to the VM system is quite challenging, 
due to the high-dimensionality (6D-|-time) of the Vlasov equation, multiple temporal and spatial scales 
associated with various physical phenomena, nonlinearity, and the conservation of physical quantities due 
to the Hamiltonian structure [43l |44] of the system. Particle-in-cell (PIC) methods [3 [35] have long been 
very popular numerical tools, in which the particles are advanced in a Lagrangian framework, while the field 
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equations are solved on a mesh. This remains an active area of research [22j . In recent years, there has 
been growing interest in computing the Vlasov equation in a deterministic framework. In the context of the 
Vlasov-Poisson system, semi-Lagrangian methods [TTl [5S] , finite volume (flux balance) methods HSl [21] , 
Fourier- Fourier spectral methods ,39i,40). and continuous finite element methods 58, 59j have been proposed, 
among many others. In the context of VM simulations, Califano et al. have used a semi-Lagrangian approach 
to compute the streaming Weibel (SW) instability [9] , current filamentation instability [42] , magnetic vortices 
[5], magnetic reconnection j7j. Also, various methods have been proposed for the relativistic VM system 

[51 H [Ml [35] . 

In this paper, we propose the use of discontinuous Galerkin (DG) methods for solving the VM system. 
What motivates us to choose DG methods, besides their many widely recognized desirable properties, is that 
they can be designed systematically to be as accurate as one wants, meanwhile with provable conservation 
of mass and possibly also the total energy. This is in general hard to achieve within other numerical method 
frameworks for simulating the VM system. The proposed scheme employs DG discretizations for both 
the Vlasov and the Maxwell equations, resulting in a consistent description of the distribution function 
and electromagnetic fields. We will show that up to some boundary effects, depending on the size of the 
computational domain, the total charge (mass) is conserved and the total energy can be preserved with a 
suitable choice of the numerical flux for the Maxwell equations and underlying approximation spaces. Error 
estimates are further established for several flux choices. The DG scheme can be implemented on both 
structured and unstructured meshes with provable accuracy and stability for many linear and nonlinear 
problems, it is advantageous in long time wave-like simulations because it has low dispersive and dissipative 
errors W, and it is very suitable for adaptive and parallel implementations. The original DG method was 
introduced by Reed and Hill |51j for a neutron transport equation. Lesaint and Raviart j41j performed 
the first error estimates for the original DG method, while Cockburn and Shu in a series of papers [HI [III 
rroi \W[ [inj developed the Runge-Kutta DG (RKDG) methods for hyperbolic equations. RKDG methods 
have been used to simulate the Vlasov-Poisson system in plasmas [Ml [331 [H] and for a gravitational infinite 
homogeneous stellar system [T^. Some theoretical aspects about stability, accuracy and conservation of 
these methods in their semi-discrete form are discussed in [33l[3l[2]. Recently, semi-Lagrangian DG methods 
[521 [50] were proposed for the Vlasov-Poisson system. In [37l[38], DG discretizations for Maxwell's equations 
were coupled with PIC methods to solve the VM system. 

The rest of the paper is organized as follows: in Section [2| we describe the numerical algorithm. In 
Section |3] conservation and the stability are established for the method. In Section |4] we provide the error 
estimates of the scheme. Section [5] is devoted to discussion of simulation results. We conclude with a few 
remarks in Section [H 

2. Numerical Methods. In this section, we will introduce the DG algorithm for the VM system. We 
consider an infinite, homogeneous plasma, where all boundary conditions in x are periodic, and /(x, ^,i) is 
assumed to be compactly supported in ^. This assumption is consistent with the fact that the solution of 
the VM system is expected to decay at infinity in ^-space, preserving integrability and its kinetic energy. 

Without loss of generality, we assume il^ = {—L^, L^Y^ and = [— i^, L^Y^ , where the velocity space 
domain Q,^ is chosen large enough so that / = at and near the phase space boundaries. We take = = 3 
in the following sections, although the method and its analysis can be extended directly to the cases when 
dx and take any values from {1,2,3}. 

In our analysis, the assumption that /(x,^,t) remain compactly support in ^, given that it is initially 
so, is an open question in the general setting. The answer to this question is important for proving the 
existence of a globally defined classical solution, and its failure could indicate the formation of shock-like 
solutions of the VM system. Whether or not the three-dimensional VM system is globally well-posed as a 
Cauchy problem is a major open problem. The limited results of global existence without uniqueness of weak 
solutions and well-posedness and regularity of solutions assuming either some symmetry or near neutrality 
constitute the present extent of knowledge [2^ 1301 [^ 1^ E51 E5ll?7] . 

2.1. Notations. Throughout the paper, standard notations will be used for the Sobolev spaces. Given 
a bounded domain _D g M* (with ★ = dx,d^, or dx + d^) and any nonnegative integer m, H"^{D) denotes 
the L^-Sobolev space of order m with the standard Sobolev norm || • and W"^''°°{D) denotes the L°°- 

Sobolev space of order m with the standard Sobolev norm || • ||m,oo,_D and the semi-norm | • |„i,oo,_D- When 
m^O, we also use H^{D) = L^{D) and W^'^^iD) = L^{D). 
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Let Ti^ — {Kx} and = {K^} be partitions of fl^ and fl^, respectively, with and being (rotated) 
Cartesian elements or simpliccs; then Th = {K : K = x K^,\/Kx £ Ti^,\/K^ G 7^f} defines a partition 
of il. Let £x be the set of the edges of T^f and £^ be the set of the edges of T^; then the edges of Th will 
be £ = {Kx X : VK^ G T^.^e^^ G £^} U {ex x iiTj : Vcx G £x,^K^ G 7^}. Here we take into account 
the periodic boundary condition in the x-direction when defining £x and £. Furthermore, £^ = £| U £^ 

with ^1 and £| being the set of interior and boundary edges of T^, respectively. In addition, we denote 
the mesh size of 7^ as /i = max(/ia:, /i^) = maxxeTd where hx = TuayiK^^Ti^ with hx^ = diam(iir2:), 
= max^^g^c hx;^ with = diam(if^), and Hk = majc^hx^, hpc^) for K = i^^; x K^. When the mesh is 

refined, we assume both ^ and ^ '''^ are uniformly bounded from above by a positive constant do- Here 



hx.min — niinA'^.eT;== and /ij.min = ™i^i<-jg7-« ^ifj- It is further assumed that {7^*}/i is shape-regular with 
* = X or ^. That is, if px^ denotes the diameter of the largest sphere included in K^,, there is 

^ < a„ VA', G n 

for a positive constant cr* independent of h^,. 
Next we define the discrete spaces 

Gl = {5 G L^{n) : g\K=K^^K, e PHKx X K^),^Kx G T^,^K^ € , (2.1a) 
- {.g G L2(r!) : g\K G P\K),'iK G T^J , 

= {U G [L2(f],)]'^== : U|k. e [P'-(if.)]'^%Vi^, G T,:) , (2.1b) 

where P^{D) denotes the set of polynomials of total degree at most r on D, and k and r are nonnegative 
integers. Note the space Gfi, which we use to approximate /, is called P-type, and it can be replaced by the 
tensor product of P-type spaces in x and f , 

[g G : g\K=K^^K, E P^Kx) x P\K^),^Kx G Tf^^K^ e T^} , (2.2) 

or by the tensor product space in each variable, which is called Q-type 

{g G L\n) : g\K=K^^K, e Q^iKx) x Q'^{K^).yKx G r^^Vif^ G T,^} . (2.3) 

Here Q'^{D) denotes the set of polynomials of degree at most r in each variable on D. The numerical methods 
formulated in this paper, as well as the conservation, stability, and error estimates, hold when any of the 
spaces above is used to approximate /. In our simulations of Section [s] we use the P-type of (2.1a I as it is 



the smallest and therefore renders the most cost efhcient algorithm. In fact, the ratios of these three spaces 
defined in ^J^, ^ and ^ are ELo ("m-T") ■ (ELo ■ + l)"" ^it^ dx = = d. 



For piecewise functions defined with respect to 7^^ or , we further introduce the jumps and averages 
as follows. For any edge e = {Kx ^ ^x} ^ ^x, with as the outward unit normal to dK^ , g"^ = ffl^-t, 
and = U|^±, the jumps across e are defined as 

= g+n++5-n-, [U], = U+ • n+ + • n", [U], = U+ x n+ + x 

and the averages are 

{9h = l{g^+g-), {U}. = ^(u+ + u-). 

By replacing the subscript x with ^, one can define [g]^, [U]^, {g}^, and {U}^ for an interior edge oiT^ 
in £^. For a boundary edge e £ £^ with being the outward unit normal, we use 

[g]^ = gn^, {.?}« = ^5, {U}^ = . (2.4) 



This is consistent with the fact that the exact solution / is compactly supported in ^. 
For convenience, we introduce some shorthand notations. 



where again ★ is x or ^ In addition, ||5||o,£ = (llffllg ^^^^^e + WgWl^r^^xs^V^^ "^'^^^ 



\ 1/2 / 1/2 

and ||5||o,£r = ( g^dsx) ■ There are several equalities that will be used later, which can be easily 
verified using the definitions of averages and jumps. 

\{9% = {g}*[9]*^ with ★ = X or ^ , (2.5a) 
[U X V], + {V}, • [U], - {U}, • [V], = , (2.5b) 

[U X V]^; + V+ • [U]^ - U- • [V]^ = 0, [U X V]^ + V- • [U]^ -U+ • [V]^ = . (2.5c) 

We end this subsection by summarizing some standard approximation properties of the above discrete 
spaces, as well as some inverse inequalities [M]. For any nonnegative integer m, let 11™ be the projection 
onto CJ™, and 11™ be the projection onto U™, then 

Lemma 2.1 (Approximation properties). There exists a constant C > 0, such that for any g G iJ™+^(r2) 
and U e [i?™+i(f]^)]'^-, the following hold: 

\\g - n™g||o,K + h]l'\\g - U"^g\\o,aK < Ch'i^+'\\g\\m+i,K , ^KeTn, 

||u - n™u||o^K, + ||u - n™u||o.aK. < c/^™+1||u||,„+i,k., vk, e T,: , 
||U-n™u||o,oo,K. <C/i™+'||U|U+i,o,,K., VK^eTff , 

where the constant C is independent of the mesh sizes Hk and Hk^ , hut depends on m and the shape regularity 
parameters Ux and of the mesh. 

Lemma 2.2 (Inverse inequality). There exists a constant C > 0, such that for any g G P"^{K) or 
P™(if:r) X P™(is:^) with K ={KxX K^) e Th, and for any U G [P"'{Kx)]''- , the following hold: 



|V,,5l|o,K < Chj,l\\g\\o,K, llVfsllo.K < Chj,l\\g\\o,K, 



||U||o,oo,K. < Ch/;^'\\V\\o^K., mkaK. < C^V'/'l|U||o,K. , 

where the constant C is independent of the mesh sizes Hk^ , hx^ , but depends on m and the shape regularity 
parameters and of the mesh. 



2.2. The Semi-Discrete DG Methods. On the PDF level, the two equations in (1.1c I involving the 



divergence of the magnetic and electric fields can be derived from the remaining part of the VM system; 



therefore, the numerical methods proposed in this section are formulated for the VM system without (1.1c I. 
We want to stress that even though in principle the initial satisfaction of these constraints is sufficient for 
their satisfaction for all time, in certain circumstance one may need to consider explicitly such divergence 
conditions in order to produce physically relevant numerical simulations |46) . 
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Given k,r > 0, the semi-discrete DG methods for the VM system arc defined by the following procedure: 
for any K ^ x £ %, look for fh G Q^, E,i, B,i G Z^^, such that for any 5 G tj^', U, V G Z^^, 



/ dthgd^d^- [ M-V^gd^di- [ fh{Eh+(xBh)-V^gd^d^ 
JK Jk Jk 



fhS. ■ n^gdsxdS, 



dt^h ■ Udx : 



if. 



B,, • V X Udx 

A' 



J OK, 



ifhCEh +£,x Bh) ■ n^)gds^dx = , 



n,, X Bf, ■ Vds^r 



Jj, • Udx 



A' 



dtBh ■ Vrfx = - / E,, • V X Vrfx - / n, X Eft • Yds, , 



OA, 



with 



Jft(x,t)= / /ft(x,C,i)^dC 



(2.6a) 
(2.6b) 
(2.6c) 

(2.7) 



Here and are outward unit normals of dK, and SivTj, respectively. All 'hat' functions are numerical 
fluxes that are determined by upwinding, i.e., 



fh£. ■ T^x ■■ ^ M ■ ^ y{M}x + 
/ft (Eft +TxBft) • : = /ft(Eft'+T'x Bft) • 



-[fh]x ■ T^x 



(2.8a) 



{/.(Eft + e X Bft)}, + l(E/.+exBft).n,| ^^^^^^ ^ ^ ^ ^^^^^^ 



X Eft 
tlx X Bft 



X Eft = X ( {Eftj^^ + ^[Bft]^ 



n^; X Bft = ria; X ( {Bftj^; - ^[^h]T ] , 



(2.8c) 
(2.8d) 



where these relations define the meaning of 'tilde'. 

For the Maxwell part, we also consider two other numerical fluxes: central flux and alternating flux, 
which are defined by 



Central flux: Eft = {Eft}, Bft = {Bft} , 

Alternating flux: E^ = E+, B^ = B^, or E^ = E~, B^ = B^ 



(2.9a) 
(2.9b) 



Upon summing up (2.6a) with respect to K £ Th and similarly summing (2.6b I and (2.6c I with respect 
to Kx G TJf , the numerical method becomes the following: look for /ft G 5ft, Eft, Bft G UJ^, such that 



ah{fh,'E,h,'Bh; g) = , 
&ft(Eft,Bft;U,V) =;ft(Jft;U) , 

for any 5 G 5ft , U, V G Z^ft, where 

aft(/ft,Eft,Bft;g) =aft.i(/ft;5) + aft,2(//i,Eft,Bft;g) , /ft(Jft;U) = - / Jft • Udx 

6ft(Eft,Bft;U,V) = / atEft-Udx- / Bft • V x Udx - / B,, ■ [V]rds, 
Jt^ Jt^ Js^ 

dtBh ■ Vrfx + / Eft • V X Vdx + / eI • [V]rds, , 
5 



(2.10a) 
(2.10b) 



and 



ah,2ifh,'Eh,Bh;g)^- I h{Eh + ^ x Bh) ■ V^gdxdS, + f f //.(E/, + ^ x B/,) • [gj^ds^dx 



Note, ail is linear with respect to fh and g, yet it is in general nonlinear with respect to E/j and due to 



(2.8b). Recall, the exact solution / has compact support in ^; therefore, the numerical fluxes of (2.8a)-(2.8d) 



and (2.9a) and (2.9b) are consistent and, consequently, so is the proposed numerical method. That is, the 
exact solution (/, E, B) satisfies 

a,,(/,E,B;g) = 0, Vg G 
6fc(E,B;U,V) = /„(J;U), VU,VeW^^. 

2.3. Temporal Discretizations. We use total variation diminishing (TVD) high-order Runge-Kutta 
methods to solve the method of lines ODE resulting from the semi-discrete DG scheme, -^Gh = R{Gh). Such 
time stepping methods are convex combinations of the Euler forward time discretization. The commonly 
used third-order TVD Runge-Kutta method is given by 

G^^^ ^Gl + AtR{Gl) 

= \gI + ^Gf + ^Ati?(Gf ), (2.11) 

where G][ represents a numerical approximation of the solution at discrete time i„. A detailed description 
of the TVD Runge-Kutta method can be found in [S3]; see also [31] and for strong-stability-perserving 
methods. 

3. Conservation and Stability. In this section, we will establish conservation and stability properties 
of the semi-discrete DG methods. In particular, we prove that subject to boundary effects, the total charge 
(mass) is always conserved. As for the total energy of the system, conservation depends on the choice of 
numerical fluxes for the Maxwell equations. We also show that is LP' stable, which facilitates the error 
analysis of Section |4] 

Lemma 3.1 (Mass conservation). The numerical solution fh G with fc > satisfies 



d 
dt 



fhd^d^ + OhAt) ^ , (3.1) 
Th 



whe 



Qh.i{t) ^ / fhraa.x{{Eh+^xBh.)-n^,0)ds^dx . 
Equivalently, with ph{x,t) ~ J^^ fh{'x.,S,,t)d^, for any T > 0, the following holds: 

[ ph{x,T)dx+ [ QhAt)dt= ( ph{x,0)d^ . 
Jt^ Jo Jt,^ 



(3.2) 



Proof. Let (7(x,^) = 1. Note that g G Q^, for any fc > 0, is continuous and Vx5 = 0. Taking this g as 
the test function in ( |2.10a ), one has 



d 
di 



[ fhd^d^ + [ [ fh{E^+tx B,,) ■ [g^ds^d^ = 
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With the numerical flux of (2.8b) and the average and jump across f| of (2.4 1, the second term above 
becomes 



fh{'Eh+^xBh)-n^ds^dic (3.3) 
:/ / ^{{-Eh + ^xBh)-n^ + \(Eh+^xBh)-n^\)ds^d^ = Qh.i{t) , (3.4) 



and this gives (3.1). Integrating in time from to T gives (3.2). □ 

Lemma 3.2 (Energy conservation 1). For k > 2, r > 0, the numerical solution fh € 5^, 'Eh.Bh G UJ^ 
with the upwind numerical fluxes (2.8a)-(2.8d) satisfies 



d_ 
dt 



Th 



MCl'dxde + / (|E,,p + |B,,ndx + QhAt) + ^hAt) = , 



(3.5) 



with 



QhAt)= f (|[E^]r|' + |[BJ,|2)ds, , e,,,3(i)= / / MCpmax((E;,+ex B,,) •n4,0)ds^dx . 

Proof. Step 1: Let ^(x,^) ~ Note that g E Qfi for k > 2 and it is continuous. In addition, 

Vx5 = 0, V^g = 2^, and ^ x U • W^g = for any function U. Taking this g as the test function in (2.10a), 
one has 



d 
dt 



Th 



MCpdxd^-2/ hE^-^dxd^- / h{Eh + ^xBh)-m\ds^dx 



Md^ dx 
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1(E, + e X B,)/, + An, ) . i\en,)ds,d. 



2 I • J,,rfx - / / ^ ((E,, + e X B^) • + |(E,, + i x B,,) • n,|) \£,\''ds^d^ 

2/ Eh-Jhd^- [ [ fh\^\^mixx{{Eh + ixBh)-n^,0)ds^dx 
JTi' -It?" -isi 



Step 2: With U = E/^ and V = B^, \2.1Qh\ becomes 
1 d 



2dt 
1 d 
2dt 
1 d 
2(ft 
1 d 
2dt 



E,,|2(fx- / Bft-VxE,,dx-/ B,,-[Eft]^ds^ 



B;,|2dx+ / E,, •VxBftdx+/ E,,-[B^]^ds^ 



T,^ 



(|E,,|2 + |B,,|2) - / (^[Eft. X Bft], + B,, • [E,,], - E,, • [B^J, j ds^ , 

f (|E;,|2 + |B;,|2)dx+i / {\[Eh]r\' + \[Bh]r\')ds, . 



The last equality uses the formulas of the upwind fluxes (2.8c)-(2.8d) as well as (2.5b). 
Combining the results in previous two steps, one concludes (3.5). □ 

Corollary 3.3 (Energy conservation 2). For fc > 2, r > and the numerical solution fh £ G^, 
E/j,B;i G UJ^ with the upwind numerical flux (2.8a)-(2.8b) for the Vlasov part, and with either the central or 
alternating flux of (2.9a) -(2.9b) for the Maxwell part, the following holds: 



d_ 
dt 



^ fh\S.?d^di + ^^(|E/.|' + |B;,ndx^ + e„,3(0 = 



Proof. The proof proceeds the same way as for Lenima[3.2[ The only difference is that here the equahties 



( 2.5b )-( 2.5c) give 



which holds for E/j and B/^ defined in the central or alternating flux in the Maxwell solver. □ 

With either the central or alternating flux for the Maxwell solver, the energy does not change due to 
the tangential jump of the magnetic and electric fields as in Lemma |3.2[ This, on the other hand, may have 
some effect on the accuracy of the methods (See Sections |4] and [5] and also [T]). 

Remark 3.4. We note that in the conservation results of Lemmas 3.1-3.2 and Corollary 3.3, the 
conservation error terms Qh.2 ^ 0, in addition, Qh,i- o.nd Qh.a depend on the numerical solution fh on the 
outflow portion of the computational boundary in ^-space, which is determined by the numerical electric and 
magnetic fields. It is easily seen that the terms Qh,i ~ 0, for i — 1,3, by choosing the computational domain 
in ^-.space to be large enough. 

Remark 3.5. Energy conservation holds as long as S Gfi. Indeed, for k < 2, the energy conservation 
results of Lemma 3.2 and Corollary 3.3 can be obtained if one replaces with gl=gl® = {5 + 

c\i\\ V^G^^VceK}. 

Finally, we can obtain the L^-stability result for fh, a result that is independent of choice of numerical 
flux in the Maxwell solver. This result will be used in the error analysis of Section [4j 
Lemma 3.6 (L^-stability of fh)- For k>0, the numerical solution fh G satisfies 



d_ 
~dt 



\fh?d^di 



Th 



\£,-^x\\[fh]x?ds^d£, 



(Eft 



-exBft)-n5l|[/ft]^|2ds^dx = 



(3.6) 



Proof. Taking g = f^ in (2. 10a I, one gets 



with 



Observe 



2dt 



fh£, ■ S/:,^fhdxdS, 



Th 



Th 



\fhfd^di \ +i?i + i?2 = , 



fhi ■ [fhlxds^dS,, i?2 = a/i^2(/h,Eft,Bft;//i) . 



Ri 



'h K^eT,^ 



Kx 



2 

fl 



dxdS^ 
dsxdS, + 



J f 



' h o. 



Ti . 



Mf'h\x + {fh^}x ■ [fh 



'T^ JF 
'h 

fh£. ■ [fhlxds^dS, , 
1 



fh^ ■ [fh]xdsxd^ , 
fhi ■ [fhlxds^d^ , 



IC • T^x\[fh]x ■ [fh]x ds^d^ , 



J F 



{-\{!l\x + {fh}x{fh\x)-i^\\i-^x\\{fh\x? ) ds^di , 



\i ■ ^x\\\f}Ax'^dsxdi 



(3.7) 



where the fourth equahty uses the definition of the numerical flux (2.8a I and the last one is due to (2.5a I. 
Similarly, 



i?2 



E 



(Eh + e X B,0 • 



--[(E„ + ^xB,)/2]5 



{A(E, + exB^)}^.[M^ 

(-;^[/^k + {^}^•[^k)•(E,.+exB;,)^ ' 



1, 



(E,, 



'E, 



2 

ex B,, 



Eh+ixBh)-n^\[h]r[h]i]ds^dx 



■n^\\[h]i?]ds^dx 



exB,0-n5l|[Md'rfs^dx, 



where the second equality is due to • (E^j + e x B/j) = and the definition of the numerical flux in (2.8b I, 
and the third equality uses (2.5a I and E/i + f x B/i being continuous in ^. With (3.7), we conclude L? 
stabflity □ 

4. Error Estimates. In this section, we establish error estimates at any given time T > for our 
semi-discrete DG methods described in Section |2.2| It is assumed that the discrete spaces have the same 
degree, i.e., fc = r, and that the exact solution satisfies / g C'^{[G,T]]H^+'^{VL) n W'^-'°^[n)) and E, B G 
C°([0,T]; [H^^'^ {^x)Y'' n \W^'°°{^x)Y'')- Also, periodic boundary conditions in x and compact support for 
/ in e are assumed. To prevent the proliferation of constants, we use A < _B to represent the inequality 
A < (constant)i3, where the positive constant is independent of the mesh size h, h^, and h^, but it can 
depend on the polynomial degree k, mesh parameters ctq, cr^ and a^, and domain parameters and L^. 



Defining Ch = f - f and e,, = IV^f-fh, it follows that f-fu = eh-Ch- Analogously, if Ch = n^E 

B -nkr, r, _E n-feTn Tn „ J _B -nk-o n iU„ "ci in -E ^E _ j -o td _i3 /-B 



E, 



= n^^B -B,eji= n^E - E,, and ej^ = H^B - B,„ then E - E,, = - and B - B,, = - C,, 
With the approximation results of Lemma |2.1[ we have 



liaiio,o</i'-+^ 



fe+i,a, 



(4.1) 



therefore, we only need to estimate Sh, and . The remainder of this section is organized as follows: 
we first state Lemmas 4T and 42 with which the main error estimate is established in Theorem |4.3|for the 
proposed semi-discrete DG method with the upwind numerical fluxes. Then, the proofs of Lemmas |4 . 1 1 and 
|4.2|will be given in subsections |4.1|and|4.2[ Lastly, for the proposed method using the central or alternating 



flux of (2.9a|-(2.9b) for the Maxwell solver, error estimates are given in Theorem 4.6 



Lemma 4.1 (Estimate of e^). Based on the semi-discrete DG discretization for the Vlasov equation of 

{\iEh+^xBh)-n^\)\\[ehk\''ds^dx 



(|2.10a|) with the upwind flux (|2.8a|-(|2.8b|), we have 
\eh\'^dxd£, 



d_ 
dt 



le • nj;||[;r,,]j;| ds^di - 



fe+i,n(l|ef Ilo,oo,f2, + ||ef ||o,oo,n^) + |/|i,oo,n(||ef ||o,n, + Ikf llo,fi,)) l|e/i||o,n 



k+i,n 



|^B||l/2 



IB 




,1/2 
1/2 



IE 



,1/2 

l0,OO,Oa 



e X B/j) • n^\\[eh]\'^ds^dx 



k+i,n 



1/2 



'dsxd^ 



(4.2) 



with 



fc+1,0 ■ 



(l + ||E||i,ooA 
(||E||fc+i,n, -Hi 



B 



||B||i_oo,n,) ||/||fc+i,a 
|fe+i,n^) |/|i,oo,o ■ 
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Lemma 4.2 (Estimate of and £^)- Based on the semi-discrete DG discretization for the Maxwell 
equations of (2.10b) with the upwind flux (2.8c)-(|2.8d|), we have 



d_ 
dt 



[ (|sfp + |sfp)rfx 
<{\\eh\Wn + h'^+^ 



(|[£f]r| 



(4.3) 

1/2 



Theorem 4.3 (Error estimate 1). For k >2, the semi-discrete DG method of ( |2.10a[ )- p.l0b[ ), for the 

(4.4) 



Vlasov-Maxwell equations with the upwind fluxes of (2.8a)- p.8d ), has the following error estimate 
\\{f~h){t)\\ln 



(E - ^h)mi,n. + ll(B - ^h)mt,n. < Gh 



2k+l 



V t e [0, T] . 

Here the constant C depends on the upper hounds of \\dtf\\k+i, a, ||/|U+i,n, |/|i,oo,n, ||E||i,oo,n,, ||B||i^oo,n^; 
||E||fc+i^sij;, I |B| Ifc+i^Q^ over the time interval [0,T], and it also depends on the polynomial degree k, mesh 
parameters ao,ax and a^, and domain parameters and L^. 

Proof. With several applications of Cauchy-Schwarz inequality and 



fe+l,f2 



1 



lEl 



,1/2 

l0,oo,O^ 



IBI 



,1/2 



Eq. (|4.2|) becomes 

\eh\'^dxd£, 



d_ 
dt 



Th 



<c (h'^+^K^ + (/i^||/|U+i,o(||£f ||o,oo,n. + Ikf llo,co,oj + |/|i,oo,a(||£f ||o,o. + Ikf llo,oj)^ 



fc+l,o(lkf llo,oo,0^ 



)) + \\eH\\ln 



<c [h''+'A' + h^'il + /»)||/||ti.o(ll^f llo.oo,o. + Ikf llo,oo.nj + I/I?. 
+ \\eH\\h 



ll^^f llo.n^)) 



Here and below, the constant c > only depends on fc, mesh parameters crQ,a.j. and a^, and domain 
parameters and L^. Moreover, with the inverse inequality of Lemma 2.2 ^nri — ^ 
bounded by ctq when the mesh is refined, we have 



and 



'^^''(lls^f llo,oo,a^ + Ikf llo,oo,a^ 



< ch 



2k- 



l^f llo,o^ 



being uniformly 
(4.5) 



and this leads to 



d 
dt 



Th 

<c ( /i^^+ip + [h^'^-^-^il + h)\\f\\l^,,, + |/|?.^.J(||£f IIL.. + Ikf ll^,o J ) + WenWin 



(4.6) 



Recall dj. = 3, then for fc > 2, there \s 2k — d^ > and therefore /i^*^ '^^ < oo. Similarly, with the 

(4.7) 



Cauchy-Schwarz inequality, (4.3) becomes 



d^ 
dt 



<c{\\e,\\l^ + h''^+^^ 



k+i,n 



h: 



2k+l/ 



|E||fc+i,o, 



Now, summing up (4.6) and (4.7), we get 



d 
dt 



\ehVd^di + 



7? 
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e 



|B||fc+i_o,)) 



l^f llo.n^ 



Th 



-S|2 



c?x 



Here A depends on (/,E,B) in their Sobolev norms \\dtf\\k+i,n, ||/|U+i,o, |/|i,oo,n, ||E||i,oo,o,, | |B| |i_oo,o, , 
||E||fc+i.n^, ||B||fe+i_n^ at time t, and 6 depends on ||/||fc+i,n and |/|i,oo,n at time t. Both A and Q depend 
on the polynomial degree k, mesh parameters ao,ax and a^, and domain parameters and i^. Now with 
a standard application of Gronwall's inequality, a triangle inequality, and the approximation results of (4.1 1, 
we conclude the error estimate (4.4). □ 



Remark 4.4. The 



n 4-3 shows that the proposed methods are (fc + ■^)-th order accurate, which is 
standard for upwind DG methods applied to hyperbolic problems on general meshes. The assumption on the 
polynomial degree fc > 2 is due to the lack of the L°° error estimate for the DG solutions to the Maxwell 



solver and the use of an inverse inequality in handling the nonlinear coupling (see ( |4.5| -(4.7) in the proof of 
Theorem 4-3). If the computational domain in x is one- or two-dimensional (d^ = 1 or 2), then Theorem 



4-3 holds for fc > 1. 



If the upwind numerical flux for the Maxwell solver ( 2.10b I is replace d by either the central or alternating 
flux (2.9a)-(2.9b|, we will have the estimates for ef and in Lemma 4.5 instead, provided an additional 



assumption is made for the mesh when it is refined. That is, we need to assume there is a positive constant 
S < 1 such that for any Kx , 



5 < 



h 



< 



1 



(4.8) 



where K^' is any element in Tf^ satisfying K^' n ^ 0. 

Lemma 4.5 (Estimate of and with the non-upwinding flux). Based on the semi-discrete DG 



discretization for the Maxwell equations of (2.10b), with either the central or alternating flux of (2.9a 



(2.9b), we have 



k-\-i,n)\\Sh \\o,n^ 



■ c{d)hlm\k+i,n. 



(4.9) 




1/2 



£-|^)dx 



The proof of this Lemma is given in Subsection 4.3 With Lemma 4.5 and a proof similar to that of 
Theorem |4.3[ the following error estimates can be established, but the proof is omitted. 

Theorem 4.6 (Error estimate 2). For fc > 2, the semi-discrete DG method of (|2.10a|)-(|2.10b|) for 



Vlasov-Maxwell equations with the upwind numerical flux (2.8a) -(2.8b) for the Vlasov solver and either the 



central or alternating fluxes of (2.9a)-(2.9b) for the Maxwell solver, has the following error estimate: 



IK/ - h){t)\\ln + IKE - ^H)mU. + ll(B - Bn)mln. < Ch 



2k 



V i e [0, T] 



(4.10) 



Besides the dependence as in Theorem 4-3. the constant C also depends on 5 of (4 



Theorem |4.6| indicates that with either the central or alternating numerical flux for the Maxwell solver, 
the proposed method will be fc-th orde r accurate. Also, one can see easily that the accuracy can be improved 
to (fc + 5)-th order as in Theorem 4.3 if the discrete space for Maxwell solver is one degree higher than that 
for the Vlasov equation, namely, r = fc + 1. This improvement will require higher regularity for the exact 
solution E and B. 

In [2 , optimal error estimates were established for some DG methods solving the multi-dimensional 



Vlasov-Poisson problem on Cartesian meshes with tensor-structure discrete space, defined in (2.3), and 



fc > 1. Some of the techniques in [2] are used in our analysis. In the present work, we focus on the P-type 



space Qi in (|2.1a|) in the numerical section, as it renders better cost efficiency and can be used on more 



general meshes. Our analysis is established only for fc > 2 due to the lack of the L°° error estimate of the 



DG solver for the Maxwell part which is of hyperbolic nature, as pointed out in Remark 4.4 



In the next three subsections, we will provide the proofs of Lemmas 4.1 4.2 and 4.5 
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4.1. Proof of Lemma |4.1[ Since the proposed method is consistent, the error equation is related to 
the Vlasov solver, 



ah(/,E,B;5r,0 - ah{fh,'Eh,Bh; gn) = 0, ^gn G Qh 



Note, Eh, & G^' by taking gh = Sh in {A.ll), one has 

ah{£h,'Eih,Bh',£h) ~ a/i(n'^/, Eh,Bh',£h) — a-hif, E, B; Eh) ■ 
Following the same lines as in the proof of Lemma |3.6[ we get 



a.h{£h, E/i, B/i; Eh) 



1 d 
'2dt 



n 



IC • rix\\[£h]x\'^dsxd^ 



(4.11) 
(4.12) 

(4.13) 



If [ |(E„+exB^).n^||M^|2ds^dx 



Next we will estimate the remaining terms in (4.12). Note 

a,,(nV, E/,, B;,; e,,) - ah{f, E, B; Sh) = Ti + , 

where 

Ti = ah,ii'n''f;eh) - a/i,i(/;e/0 = aft,i(C/i; ^h) , 
1^2 =aft,2(nV,Eft,Bft;eft)-a/,,2(/,E,B;e;,) . 

Step 1: estimate of Ti. We start with 



Ti - / (9ta)£/.rfxde - / a^ Vxeh^xde 



J F 



ChC-[£h]xds.,d^^Tn+Ti2+T, 
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It is easy to verify that cJill'^ = H'^dt, and therefore dtCh = ^^{dtf) ~ (dtf)- With the approximation result 
of Lemma 2.1 we have 

1^11 1- / {dtCk)ehd^d^ 



n 



< \\dtCh\\o,n\\^h\\o.n < /j'+'ll<9t/||fe+i,n||e„||o,n 



(4.14) 



Next, let ^0 be the projection of the function ^ onto the piecewise constant space with respect to 7^, 



then 



Ti2 = - / a(C - • Vxe/.dxd^ - / a^o • y^ehdKd^ 



(4.15) 



Since ' Vx^/i € and C/i = n'^/ — / with II'^ being the projection onto the second term in (4.15) 
vanishes. Hence 



ri2| < / ia(^-6)-Vx£/.|dxde 



ITh 

0,00, 0,K) 

<||e-Co||o,oo,n, E /^^+'/^kJI/IU+i, 

< /icllClkoo,n5/i''||/IU+i,n||£h||o,n 
</i'+^||/lk+i,n||£.||o,n . 



(4.16) 



The third inequality above uses the approximating result of Lemma 2.1 and the inverse inequality of Lemma 
2.2 The fourth inequality uses an approximation result similar to the last one of Lemma 2.1 and — — 
being uniformly bounded by ctq when the mesh is refined. 
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Next, 



Ti3 = / , / ( + ^^[CJ. ) • [shUs^d^ 



J f 



I . {Ch}x{^ ■ nx)nx + ^-^-^[Chlx ] ■ [£h]xds^d^ , 
where iix is the unit normal vector of an edge in £x with either orientation, that is xyx = n^,, or — n^;. Then, 

WChU 



|Ti3|< I j (l^-n^KKaWH- 



\[£h]x\dsxd(, 

1/2 / s 1/2 



< 



(^j^^ ^ 2(|{C4.P + ■ r^x\ds,d^ (^J^^ ^ I? • nMshU'ds^d^ 



/ \ 1/2 

<ll^llolnJiailo,r,^xf. [I^J^ \i-n.Mek]x\'dsxd^] 



\ 1/2 

<h'+'^\\f\\k+i,n{ I I \^-nx\\[eh]x\'dsxdA . (4.17) 



The approximation results of Lemma |2.1| are used for the last inequality. 

Step 2: estimate of T2. Note, 

T2 = a/.,2(n'=/, E,„ B,,; Eh) - atM^ E, B; Eh) 

= 0'h,2{Ch, E/i, Bft ; e;i) + ah^2{f, Eh,^h',£h) — ah,2{f, E, B; Sh) = T21 + T22 + 123 , 

with 

T2i = - [ ChiEh + e X Bft) • V^Ehdxdt T22 = I I C/i(Eh'+rx B,,) • [Eh]^ds^dx, 

723 — a/i,2(/, E/i, B/i; Eh) — ah,2{f, E, B; Eh) ■ 

For T21, we proceed as for the estimate of T12. Let Eq = n°E, Bp — n°B be the projection of E, B, 
respectively, onto the piecewise constant vector space with respect to Th, then 

/ Ch{Eh +^xBh)- V^Ehdxd^ = [ ChiEh - Eo + e X (B^ - Bq)) • V^e^dxd^ 

+ Oi(Eo + e X Bo) • V^Ehdxde , 
and the second term above vanishes due to (Eq + ^ x Bq) • V^e/i G t/^, and therefore 

I / Chi'Eh + e x B,,) • \/^Ehdxd(\ < ( |a(Eft - Eo + e X (Bft - Bq)) • V^e/Jdxd^ , 

< (||E;,-Eo+ex (B;,-Bo)||o.oo,o) {hK\\\Ch\WK){hK,\W ^Eh\\o,K) . 

K^-x.K^=KeTh 

<(||E„-Eo||o.oo.o. + ||(B„-Bo)||o,oo,nJ h)+^h-K\\\f\\k+i,K\\^h\WK , 

K^xK(=K(^Th 

< /i'||/IU+i,o(||ef ||o,oo,n. + lief ||o,oo,o. + ||n^E - Eo||o,oo,n. + ||n^B - Bo||o,oo,nJ||eh||o,n ■ 
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Note that H^E - Eq = n^(E - Eq), and U'^ is bounded in any iP-norm (1 < p < oo) [20l E], then 
and similarly ||n^B - Bo||o,oo,n, < /j2;||B||i,oo,o,- Hence 

(4.18) 

</»"||/||fc+i,o(||e^'||o.oo,a. + \\£h\koo^n^ + /ix(||E||i,oo.n. + ||B||i,oo,oJ)||£/.||o,o ■ 
For T22, we follow the estimate of Ti^. Note that E^ and only depends on x, and ^ is continuous, 

4(E,, + 5 X B^) • [eh]ids^dx\ 



[ a(E/.+exB;,)-v4e,,dxde 



{a(E. + exB.)}, + &±i4M:i^[c,] 



{a}e((E, + e X B,) . h^)h^ + l(E/.+exB,).ne| ^^^^j^^ ^ [e^J^ds^da:! , = or - 



< 



/ / (m,+^xB,)-n^\{\{a}^\ + \^^\))\[s,]^\ds^dx, 

J^J^2\(E,,+^xB,,)-n^\\{Cl}\ds^dx^ (^j^J^ {(Eh + ^ x B,,) ■ n^\\[eh]\^ds^dx 



1/2 



< 



In addition, 



/ / 2\{E^+^xB^)-n^\\{Cl}\ds^dx 



1/2 



/r- J£, 

< IITP., _L C -R, 



|Eh + ^xB;,||;/4^||a||o,r,^x£, 
^ IIC/i||o,r^-xf5(||E/i||o^^ + ||B,i||o^^ 

</^'+^||/lk.+i.f2(||£f|lolo. + lkfllolo. 



lEl 



,1/2 

l0,CX3,£l3; 



IRI|l/2 ^ 



and therefore 



r22</i'=+^||/IU+i,o(||£f|loln. + lk 



S||l/2 

/i llo.oo.ri^ 



lEl 



,1/2 



IBI 



,1/2 



/ / \{Y.i,+ixBy,)-n^\\[eh]\'ds^dx\ 
J£, J 



a,oo,n^ ^ I 10,00,0^ 
1/2 



(4.19) 



Finally, we estimate r23. Since / is continuous in ^, and • (E/j — E + ^ x {B^ — B)) = 0, 
723 — ah,2{f, E/i, B/;; £/i) — a/,.^2(/, E, B; £h) 

= - j /(E,, - E + ^ X (B,, - B)) • V^Ehd^di + I j f{Eh - E + ^ x (B^ - B)) • [ehkds^dx 

= [ Vj/ • (E,, - E + ^ X (B,, - B))e,,dxde ; 



therefore, 



IT23I < ||E,, - E + e X {Bh - B)||o,n|/|i,oo,n||e/.||o,n , 

< (||E^ - E||o,fi, + ||(B,, - B)||o,nJ|/|i,oo,o||eh||o.n , 

< (ll^f llo.n. + lief llo.n. + llCf Ilo,o. + ||Cf Ilo,oj|/|i,oo,n||£.||o,n , 

< (Ikf llo.fi. + WefAkn. + /i^^(||E|U+i,o. + ||B|U.+i,oJ)|/|i,oo,o|kh||o,o 



(4.20) 



Now we combine the estimates of (4.14) and (4.16l-(4.20l, and get the result of Lemma 4.1 
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4.2. Proof of Lemma |4.2[ Since the proposed method is consistent, the error equation is related to 
the Maxwell solver, 



&^(E-E,,B-B„;U,V) J,„U), yv,veu^ 

Taking the test functions in (4.21) to be U = ef and V = gives 



E /-B, R 



Following the same lines of Step 2 in the proof of Lemma 3.2 



L r^E ^B. ^B\ 



-B\2\ 



dx+ - 



-El |2 



dSr, 



(4.21) 



(4.22) 



(4.23) 



It remains to estimate the two terms on the right side of (4.22), 



, f^E >B. E ^B\ 



C^Vxefdx- 

Ch ' [£h\TdSx , 



< 



'7~a; / '7~a; 

' ' h 

1/2 



(4.24) 



K^l'ds, 



-El |2 
'h\r\ 



sHU'ds., 



1/2 



+ ||B|U+i,oj(^ |[£f].p + |[ef].pd.. 



1/2 



</i'+^(||E||fc+i,a 



1/2 



All of the volume integrals of \A.2A\ vanish due to dfU^ = nj:at and ef ,ef , V x ef , V x ef e ^'J. And, 
for the last two inequalities, the definition of the numerical fluxes are used together with the approximation 
results of Lemma |2.1| Finally, 



,(J-J„;£f)| = | / (J - J;,) • ef dx| , 

< ||J-J,||o,nJ|£f||o,a. = 11 / (/-A)edeilo,nJ|£f||o,n. , 



< ||/-Mlo,a||ClknJ|£flka. 



< 



i\\eH\\o,n + ||a|lo,n)||£f llo.n. < (lk/.||o,a + /i'=+^||/IU+i,a)||£f ||o,n. 



(4.25) 



Combining (4.23 )-( 4.25 ), we conclude Lemma 4.2 



4.3. Proof of Lemma |4.5[ The proof proceeds in a manner similar to that of Lemma [42| of Subsection 
4.2 Based on the error equation (4.21), related to the Maxwell solver with some specific test functions, we 



get (4.22). With either the central or alternating flux of (2.9a)-(2.9b), we have 

1 d 



L (^E ^B. ^E ^B^ 



2dt 
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The same estimate as that of (4.25 1 can be obtamed for the second term on the right of (4.22 ). To estimate 
the first one, 



dtCh ■ 



(4.26) 



< 



E 



\ 1/2 



hKli\CH{' + \C^nds. 



*-B|2 



E / hKA\[e^]r\' + \[ef]r\')d.Sx 



1/2 



\C^r)ds 



2k 



lEl 



fe+l.-ftTx 



< c{6)hUm\k+i,n. 



IBI 



fc+l.Oj 




1/2 



hKAi^jii' + isHnds 



1/2 



E||2 



-B||2 ^1 



(4.27) 

(4.28) 
(4.29) 



As befor e, all volume integrals of ( |4.26[ ) vanish due to 9(11^ = nj;9t and £f,ef ,V x ef,V x ef £ Z^^. 
In (4.27), is any element containing an edge e. To get (4.281, we use the definitions of the numerical 
fluxes, jumps, as well as the assumption (4.8) on the ratio of the neighboring mesh elements. Here c{S) is a 
positive constant depending on S. We obtain (4.29) by applying an approximation result of Lemma [2J] and 
an inverse inequality of Lemma |2.2| From all the above, we conclude Lemma |4.5| 

5. Numerical results. In this section, we perform a detailed numerical study of the proposed scheme 
in the context of the streaming Weibel (SW) instability first analyzed in [35]. The SW instability is closely 
related to the Weibel instability of [ST], but derives its free energy from transverse counter-streaming as 
opposed to temperature anisotropy. The SW instability and its Weibel counterpart have been considered 
both analytically and numerically in several papers (e.g. |48l |9l HI [Tj [47] ) - here we focus on comparison with 
the numerical results of Califano et al. in [S]. 

We consider a reduced version of the Vlasov-Maxwell equations with one spatial variable, X2, and two 
velocity variables, and ^2, The dependent variables under consideration are the distribution function 
f{x2,^i,^2,t), a 2D electric field E — {Ei{x2,t), E2{x2,t),0) and a ID magnetic field B = (0,0,53(0:2,*)), 
and the reduced Vlasov-Maxwell system is 



ft + 6/.. + {El + 6S3)/6 + iE2 - aSs)/?, = , 
dB3 dEi dEi dBs . dE2 



dx2 



dt 



dt 



-]2 



where 



Ji = 



/CO f^oo poo poo 

/ /(a:2, a, 6,^)6^6^6, j2^ / f{x2,^i,a,m2d^id^: 



The initial conditions are given by 

1 



£^1(2:2,6,6,0) = £:2(a;2, 6,6,0) = 0, -83(2:2,6,6,0) &sin(fcoa:2) , 



(5.1) 
(5.2) 

(5.3) 

(5.4) 
(5.5) 
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which for 6 = is an equihbrium state composed of counter-streaming beams propagating perpendicular to 
the direction of inhomogeneity. Following [5], we trigger the instability by taking /3 — 0.01, b — 0.001 (the 
amplitude of the initial perturbation to the magnetic field). Here, fix = [Oj^yli where Ly — 27r//co, and we 
set = [—1.2, 1.2]^. Two different sets of parameters will be considered, 

choice 1 : S = 0.5, Wo.i = vo,2 = 0.3, fco = 0.2 
choice 2 : 6 = 1/6, uo,i = 0.5, uo,2 = 0.1, fco = 0.2 . 

For comparison, these are chosen to correspond to runs of [S]. 

Accuracy test: The VM system is time reversible, and this provides a way to test the accuracy of 
our scheme. In particular, let /(x, ^, 0), E(x, 0), B(x, 0) denote the initial conditions for the VM system and 
/(x,^,T),E(x,r),B(x,T) the solution at t = T. If we choose /(x, T), E(x, T), -B(x, T) as the initial 
condition at t = 0, then at i = T we theoretically must recover /(x, — ^, 0), E(x, 0), — B(x, 0). In Tables 



5.1 5.2 and 5.3 we show the errors and orders of the numerical solutions with three flux choices for 
the Maxwell's equations: the upwind flux, the central flux, and one of the alternating fluxes E/j = E^ and 
B/i = B^. The parameters are those of choice 1, with symmetric counter-streaming. In the numerical 
simulations, the third order TVD Runge Kutta time discretization is used, with the CFL number CcR = 0.19 
for the upwind and central fluxes, and CcA = 0.12 for the alternating flux in and cases. For P^, we take 
At = 0{Ax'^^^) to ensure that the spatial and temporal accuracy is of the same order. From Tables 



5.1 



5.2 



and 5.3 we observe that the schemes with the upwind and alternating fluxes achieve optimal (fc + l)-th order 
accuracy in approximating the solution, while for odd fc, the central flux gives suboptimal approximationof 
some of the solution components. 

Table 5.1 

Upwind flux for Maxwell's equations, errors and orders. Run to T=5 and back to T = 10. 



Space 




Mesh=20=' 


Mesh=403 


Mesh=803 


error 


error 


order 


error 


order 




f 


0.18E+00 


0.50E-01 


1.82 


0.13E-01 


1.96 




0.26E-05 


0.66E-06 


2.01 


0.16E-06 


2.01 


El 


0.21E-05 


0.68E-06 


1.61 


0.19E-06 


1.81 


E2 


O.lOE-05 


0.22E-06 


2.23 


0.22E-07 


3.29 


QlMl 


f 


0.56E-01 


0.77E-02 


2.87 


O.lOE-02 


2.92 




0.23E-06 


0.26E-07 


3.12 


0.32E-08 


3.06 


El 


0.16E-06 


0.16E-07 


3.32 


0.14E-08 


3.54 


E2 


0.16E-06 


0.22E-07 


2.90 


0.15E-08 


3.91 




f 


0.12E-01 


O.lOE-02 


3.56 


0.70E-04 


3.90 


B3 


0.97E-07 


0.23E-08 


5.37 


0.12E-09 


4.34 


El 


0.19E-07 


0.27E-09 


6.16 


0.57E-11 


5.54 


E2 


0.14E-07 


0.79E-09 


4.11 


0.16E-10 


5.64 



Conservation properties: The purpose here is to validate our theoretical result about conservation 
through two numerical examples, the symmetric case and the non-symmetric case. We first use parameter 
choice 1 as in the Califano et al. |^, the symmetric case with three different fluxes for the Maxwell's 
equations. The results are illustrated in Figure |5.1| . In all the plots, we have rescaled the macroscopic 
quantities by the physical domain size. For all three fluxes, the mass (charge) is well conserved. The largest 
relative error for the charge for all three fluxes is smaller than 4 x 10~^". As for the total energy, we could 
observe relatively larger decay in the total energy from the simulation with the upwind flux compared to the 
one with the other two fluxes. This is expected from the analysis in Section [3j In fact, the largest relative 
error for the total energy is bounded by 1 x lO"'' for the upwind flux, and bounded by 1 x 10~^ for central 
and alternating fluxes. In Figure |5.2[ we study the effect of enlarging the domain in the velocity direction. 
The growth in the total mass, as time approaches T = 200 when ilj = [—1.2, 1.2]^, implies that up to this 
time, a larger domain should be used in order for the assumption, / being compactly supported in ^, to still 
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Table 5.2 

Central flux for Maxwell's equations, errors and orders. Run to T=5 and back to T = 10. 



Space 




Mesh=203 


Mesh=403 


Mesh^SO^ 


error 


error 


order 


error 


order 




f 


0.18E+00 


0.50E-01 


1.82 


0.13E-01 


1.96 


B3 


0.13E-04 


0.85E-05 


0.66 


0.50E-05 


0.75 


El 


0.19E-05 


0.13E-05 


0.51 


0.58E-06 


1.17 


E2 


0.92E-06 


0.19E-06 


2.26 


0.20E-07 


3.24 




f 


0.56E-01 


0.77E-02 


2.87 


O.lOE-02 


2.92 




0.28E-06 


0.28E-07 


3.34 


0.32E-08 


3.15 


El 


0.18E-07 


0.56E-09 


5.00 


0.88E-11 


5.99 


E2 


0.16E-06 


0.22E-07 


2.90 


0.15E-08 


3.91 


QlMl 


f 


0.12E-01 


O.lOE-02 


3.56 


0.70E-04 


3.90 




O.lOE-06 


0.44E-08 


4.57 


0.16E-09 


4.81 


El 


0.46E-07 


0.82E-10 


9.12 


0.30E-10 


1.45 


E2 


0.14E-07 


0.79E-09 


4.12 


0.16E-10 


5.65 



Table 5.3 

Alternating flux for Maxwell's equation, errors and orders. Run to T=5 and back to T = 10. 



Space 




Mesh=20^ 


Mesh=40^ 


Mesh=80^ 


error 


error 


order 


error 


order 


QlMl 


f 


0.18E+00 


0.50E-01 


1.82 


0.13E-01 


1.96 




0.29E-05 


0.78E-06 


1.90 


0.22E-06 


1.83 


El 


0.24E-06 


0.35E-07 


2.74 


0.22E-08 


3.99 


E2 


O.lOE-05 


0.22E-06 


2.23 


0.22E-07 


3.29 




f 


0.56E-01 


0.77E-02 


2.87 


O.lOE-02 


2.92 


B3 


0.28E-06 


0.22E-07 


3.70 


0.18E-08 


3.63 


El 


0.32E-07 


0.30E-09 


6.72 


O.llE-10 


4.84 


E2 


0.16E-06 


0.22E-07 


2.90 


0.15E-08 


3.91 


QlMl 


f 


0.12E-01 


O.lOE-02 


3.56 


0.70E-04 


3.90 


B3 


O.lOE-06 


0.24E-08 


5.42 


0.12E-09 


4.36 


El 


0.98E-08 


O.lOE-09 


6.60 


0.90E-12 


6.80 


E2 


0.14E-07 


0.79E-09 


4.11 


0.16E-10 


5.64 



hold. This growth in relative error is not observed when fl^ = [—1.5, 1.5]^. On the other hand, the decay 
in the total energy with the upwind fluxes is largely due to the tangential jump terms in the electric and 
magnetic field as derived in Lemma 3.2. Therefore, enlarging the domain has little effects on this. As for 
the decay in total energy with central and alternating fluxes, we can see that enlarging the domain roughly 
reduces the error by half. The other part of the error is coming from the dissipative nature of the TVD-RK 
scheme that we have used. 

The two species VM system conserves the following expression for the total linear momentum: 

^fd^dx + J E X Bdx , (5.6) 

where the first term represents the momentum in the particles while the second that of the electromagnetic 
field. In fact, this is true for the full energy-momentum and angular momentum tensors |49j . Each component 



of the spatial integrand of (5.6 1, the components of the momentum density, satisfies a conservation law, a 



result that relies on both species being dynamic and one that relies on the constraint equations (1.1c 
being satisfied. However, in this paper we have fixed the constant ion background by charge neutrality 
and, consequently, momentum is not conserved in general. This lack of conservation does not appear to be 
widely known, but it is known that the enforcement of constraints sometimes, but not always, results in 
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the loss of conservation [45]. For example, the single species Vlasov-Poisson system with a fixed constant 
ion background does indeed conserve momentum. However, for the streaming Weibel application, it is not 
difficult to show that the following component is conserved: 



£,if d£,id£,2dx2+ / E2B3dx2 



(5.7) 



while the component P2 is not. Since conservation of Pi relies on the constraint equations and since our 
computational algorithm does not enforce these constraints, conservation of Pi serves as a measure of the 
goodness of our method in maintaining the initial satisfaction of the constraints. From Figure |5.1[ we see 
that all three fiux formulations conserve Pi relatively well, but, as expected, there is a large accumulating 
error in P2, particularly for the alternating flux case. 

Similarly, for a general VM system without constraints, the following expression for the total angular 
momentum is conserved: 



L = y X X ^ / dS,dx + y X X (E X B) dx 



(5.8) 



However, because the SW application breaks symmetry, there is no relevant component of the angular 
momentum that is conserved for this problem, but for a more general application one may want to track its 
conservation. 



Comparison and interpretation: In Figure 5.3 we plot the time evolution of the kinetic, electric, 
and magnetic energies. In particular, we plot the separate components defined by Ki = | J f£,id^id£_2dx2, 
K2 = ^ J f£,2dS,id£_2dx2, El — ^ J Eldx2, and E2 — ^ J i?|da;2. Figure (a) shows the transference of kinetic 
energy from one component to the other with a deficit converted into field energy. This deficit is consistent 



with energy conservation, as evidenced by Figure 5.1 Observe the magnetic and inductive electric fields 
grow initially at a linear growth rate (comparable to that of Table I of [S]). Saturation occurs when the 
electric and magnetic energies simultaneously peak at around i = 70 in agreement with [9]; however, in 
our case we achieve equipartition at the peak, which may be due to better resolution. Here we have also 



shown the longitudinal component £'21 not shown in [5], which in Figure 5.4 is seen to grow at twice the 
growth rate. This behavior was anticipated in [TU] in the context of a two-fluid model and seen in kinetic VM 
computations of the usual Weibel instability [17] . It is due to wave coupling and a modulation of the electron 
density induced by the spatial modulation of i?| . The growth at twice the growth rate of the magnetic field 



B3 is seen in Figure |5.4[ and the density modulation, including the expected spikes, is seen in Figure |5.5| 
We have also calculated the first four Log Fourier modes of the fields Ei, E2, -B3, and these are shown in 
Figure 5.6 Here, the n-th Log Fourier mode for a function W(x,t) [M] is defined as 



logFMn{t)= log 



10 



L 



\ 



W{x, t) shi{knx) dx 



W{x, t) cos{knx) dx 



In Figures [577| and 5.8 we plot the 2D contours of / at selected locations X2 and time t, when the upwind flux 
is used in the Maxwell solver. The times chosen correspond to those for the density of Figure |5.5[ and we 
see that at late times considerable fine structure is present, which is consistent with the Log Fourier plots. 
For completeness, we also include in Figure [5^ plots of the electric and magnetic fields at the final time for 
our three fluxes. 
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Fig. 5.1. Streaming Weibel instability with parameter choice 1 as in Califano et al. JSj (6 = 0.5, I'd, 1 = 
O.2J, the symmetric case. The mesh is 100^ with piecewise quadratic polynomials. Time evolution of mass, 
momentum for the three numerical fluxes for the Maxwell's equations. 
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Fig. 5.2. Streaming Weibel instability with parameter choice 1 as in Califano et al. [9] (S = 0.5, fo,! = fo,2 = 0.3, fco = 
0.2), the symmetric case. Effects of enlarging the domain. The runs are conducted on the same mesh size, but with different 
domains in the velocity space. The small domain is = [—1.2, 1.2]'^, the large domain is fl^ = [—1.5, 1.5]'^ 




t t 

(a) Kinetic energy (b) Electric and magnetic energy 

Fig. 5.3. Streaming Weibel instability with parameter choice 1 as in {PP (S = 0.5,i)o,i = "0,2 = 0.3, fco = 0.2). The mesh 
zs 100^ with piecewise quadratic polynomials. Time evolution of kinetic, electric and magnetic energy by alternating flux for 
the Maxwell's equations. 




Fig. 5.4. Growth rates for streaming Weibel instability with parameter choice 1 as in fP'/ (5 = 0.5, vo,l = ^^0,2 = 0.3, fco = 
0.2). The mesh is 100^ with piecewise quadratic polynomials. Time evolution of kinetic, electric and magnetic energy by 
alternating flux for the Maxwell's equations. 



22 




Fig. 5.5. Plots of the computed density function for the streaming Weihel instability, with parameter choice 1 as in 0' 
(5 = 0.5, HQ, 1 = "0,2 = 0.3, feo = 0.2j, at selected time t. The mesh is 100"^ with piecewise quadratic polynomials. The upwind 
flux is applied. 
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Fig. 5.6. Streaming Weibel instability with parameter choice 1 as in (5 = 0.5, iio,! = ''0,2 = 0.3, fco = 0.2). The mesh 
IS 100^ with piecewise quadratic polynomials. The first four Log Fourier modes of E\ , E2 , B3 computed by the alternating flux 
for the Maxwell's equations. 



24 





Fig. 5.7. 2D contour plots of the computed distribution function ff^ for the streaming Weibel instability, with parameter 
choice 1 as in (5 = 0.5,do,i = ^'0,2 = 0.3, fcg = 0.2), at selected location X2 and time t. The mesh is 100^ with piecewise 
quadratic polynomials. The upwind flux is applied. 




Fig. 5.8. 2D contour plots of the computed distribution function ff^ for the streaming Weibel instability, with parameter 
choice 1 as in (5 = 0.5,do,i = ^'0,2 = 0.3, fcg = 0.2), at selected location X2 and time t. The mesh is 100^ with piecewise 
quadratic polynomials. The upwind flux is applied. 






(e) Electric field, alternating flux (f) Magnetic field, alternating flux 
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Fig. 5.9. Streaming Weibel instability with parameter choice 1 as in l9j (S = 0.5, no, i = ^0,2 = 0.3, fco = 0.2). The mesh 
is 100^ with piecewise quadratic polynomials. The electric and magnetic fields at T = 200. 



Now consider the nonsymmetric parameter set, choice 2 of the Cahfano paper. The results for this 
choice are gathered in Figures |5.10[ |5.1H |5.12[ |5.13[ |5.14[ |5.15[ and |5.16[ which parallel those of parameter 
choice 1. Again we observe relatively larger error in the total energy when the upwind flux is used to solve 
Maxwell's equations. But overall, both mass and total energy are very well preserved. Insofar as we can 
make comparison with [9], our results are in reasonable agreement. Similar energy transfers take place, but 
the equipartition of the magnetic and electric energies at the peak is not achieved. All modes saturate now 
at nearly the same values, evidently resulting from the broken symmetry. Also, at long times, contours of 
the distribution function are displayed. Here the wrapping of the distribution function as two intertwined 
distorted cylinders is observed as in [5] , although for late times there is a loss of localization. 
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Fig. 5.10. Streaming Weibel instability with parameter choice 2 as in Califano et al. 5 = l/6,tio,i = 0.5, vo, 2 = 
0.1, fco = 0.2), the non-symmetric case. The mesh is 100^ with piecewise quadratic polynomials. Time evolution of mass, total 
energy with three numerical fluxes for the Maxwell's equations. 
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(a) Kinetic energy (b) Electric and magnetic energy 

Fig. 5.11. Streaming Weibel instability with parameter choice 2 as in 0' = l/6,i)o,i = 0.5, do, 2 = 0.1, fco = 0.2). The 
mesh is 100^ with piecewise quadratic polynomials. Time evolution of kinetic, electric and magnetic energy by alternating flux 
for the Maxwell's equations. 
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Fig. 5.12. Streaming Weibel instability with parameter choice 2 as in (& = 1/6, i)o,l = 0.5, 1)0,2 = 0.1, fco = 0.2). The 
mesh is WQp with piecewise quadratic polynomials. The first four Log Fourier modes of Ei, E2, S3 computed by the alternating 
flux for the Maxwell's equations. 
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Fig. 5.13. Plots of the computed density function pf^ for the streaming Weibel instability, with parameter choice 2 as in 
fpy (5 = 1/6,1)0,1 = 0.5, "Uq, 2 = 0.1, fco = 0.2), at selected time t.The mesh is 100'' with piecewise quadratic polynomials. The 
upwind flux is applied. 




Fig. 5.14. 2D contour plots of the computed distribution Junction fi^ for the streaming Weibel instability, with parameter 
choice 2 as in (S = 1/6, no.l = 0-5, vo,2 = 0.1, fco = 0-2 ), at selected location X2 and time t. The mesh is 100^ with piecewise 
quadratic polynomials. The upwind flux is applied. 




Fig. 5.15. 2D contour plots of the computed distribution Junction fi^ for the streaming Weibel instability, with parameter 
choice 2 as in (S = 1/6, no.l = 0-5, vo,2 = 0.1, fco = 0-2 ), at selected location X2 and time t. The mesh is 100^ with piecewise 
quadratic polynomials. The upwind flux is applied. 




(a) Electric field, upwind flux (b) Magnetic field, upwind flux 




(c) Electric field, central flux (d) Magnetic field, central flux 




(e) Electric field, alternating fiux (f) Magnetic field, alternating flux 



Fig. 5.16. Streaming Weibel instability with parameter 3wice 2 as in (3' = 1/6, rio.l = 0.5,tio,2 = 0.1, fco = 0.2). The 
mesh is 100'^ with piecewise quadratic polynomials. The electric and magnetic fields at T = 200. 



6. Concluding Remarks. In summary, we have developed discontinuous Galerkin methods for solving 
the Vlasov-Maxwell system. We have proven that the method is arbitrarily accurate, conserves charge, can 
conserve energy, and is stable. Error estimates were established for several flux choices. The scheme was 
tested on the streaming Weibel instability, where the order of accuracy and conservation properties were 
verified. 

In the future, we will explore other time stepping methods to improve the efficiency of the overall 



algorithm. In our development, the constraint equations of (1.1c) were not considered; in the future, we 
plan to investigate them together with some correction techniques for the continuity equation. The proposed 
method has been clearly established as sufficient for investigating the streaming Weibel instability, and the 
long time nonlinear physics of this system can be further investigated and modeled. In the future, we will also 
apply the method to study other important plasma physics problems, especially those of higher dimension. 
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